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We introduce the Poisson-Boltzmann cell model for spherical colloidal particles with a hetero- 
geneous surface charge distribution. This model is obtained by generalizing existing cell models 
for mixtures of homogeneously charged colloidal spheres. Our new model has similar features as 
Onsager's second-virial theory for liquid crystals, but it predicts no orientational ordering if there 
is no positional ordering. This implies that all phases of heterogeneously charged colloids that are 
liquid-like with respect to translational degrees of freedom are also isotropic with respect to particle 
orientation. 



I. INTRODUCTION 

Already a long time ago Marcus [l|, Ohtsuki et al. 
and Alexander et al. [3| realized that the cell model 
approach of Wigner and Seitz [3| to calculate the prop- 
erties of electrons in solids can also be applied to col- 
loidal matter. In this case there are no quantum effects, 
and in stead of a wave function one calculates the ion- 
distributions around charged colloidal particles dispersed 
in a liquid medium. In the simplest case, the colloidal dis- 
persion consists of one colloidal species immersed in a 1:1 
electrolyte solution. The colloidal particles are homoge- 
neously charged and have a spherical shape. Addition- 
ally, the simplification of taking a spherical Wigner-Seitz 
cell — instead of space filling — is justified in fluid phases 
with no broken translational symmetry. 

A number of extensions has been made to this basic 
cell model. There is the eccentric Poisson-Boltzmann 
cell model 0, and the heterogeneous — or polydis- 
perse — cell model [f| 0|, to describe mixtures. Addi- 
tionally, the cell model has been extended by applying 
cylindrical Wigner-Seitz cells for the description of disc- 
shaped and rodlike particles d,Q. In the present paper, 
we make an extension towards particles with an heteroge- 
neous surface charge distribution, such as patterned col- 
loids [l0( or Janus particles [ll| . Janus particles are char- 
acterized by two distinct regions of surface area. Each of 
these "faces" has a different chemical composition, which 
can create spontaneous aggregation [T3|. Our aim is to- 
wards a description of self-organization of these particles 
from single-particle properties. The cell model can be a 
powerful tool, giving a simple description of these com- 
plex systems. From this description, a number of ther- 
modynamic quantities can be derived such as free energy, 
and osmotic pressure. 

We investigate how surface charge heterogeneity in- 
fluences the distribution of particle orientations in the 
case of a homogeneous positional distribution. This case 
can be considered as a simplified description of a fluid of 
these particles, but also as an over-simplified description 
of a solid. Since the particle interactions are implemented 
through the boundary conditions of the Wigner-Seitz cell, 
one can choose to neglect the correlations between the 
orientational and positional ordering. The basis of our 



model is the generalization of the Poisson-Boltzmann cell 
model by Biesheuvel et al. Q considering mixtures, to- 
gether with the insight of Onsager [l3| that a distribution 
of orientations can be considered in the same way as mix- 
tures. In principle, this model can be used to predict the 
phase behavior of a large class of colloidal and nanoparti- 
cles, because an anisotropic (as well as spherical) particle 
shape can be treated. However, here we restrict ourselves 
to a system of spherical colloidal particles. 




Figure 1: Illustration of the Poisson-Boltzmann cell model. 
A colloidal particle of radius a is surrounded by ions. It is 
situated in the center of a spherical cell of radius R. 



II. THE CELL MODEL FOR MIXTURES 

We start by giving the prerequisites for the ordinary 
spherically symmetric cell model (see Fig.[T]). In this sim- 
ple case, we consider a system of N identical, spherically 
symmetric colloidal particles in a volume V . These par- 
ticles have radius a and a surface charge Ze, where e is 
the elementary charge. The system is presumed to be in 
osmotic contact with a reservoir of monovalent cations 
and anions at total concentration 2p s , with charge ±e. 
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Both the particle and the surrounding medium are con- 
sidered to be dielectric media, where e ln and e out are the 
dielectric constants of the particle and the solvent, re- 
spectively. 

Within Poisson-Boltzmann theory, one relates the 
denisty profiles p± (r) of the ions to the electrostatic po- 
tential \&(r) in a fixed configuration of the colloidal parti- 
cles. This complicated nonlinear iV-body problem can be 
simplified considerably by considering a single colloidal 
particle in the center of a Wigner-Seitz cell, surrounded 
by the ions. This cell is assumed to have a spherical 
shape (of radius R) instead of space filling. The volume 
of the cell is fixed by the available volume per particle, 
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We treat the ions in a mean field description, such that 
we obtain the Poisson-Boltzmann (PB) equation 



charge Zie. Again, the particles have radius a, and are 
positioned in the center of a spherical Wigner-Seitz cell 
of radius R. The PB-equation (|2|) is solved separately for 
each colloidal species i, and the notion of electroneutral- 
ity is applied to the system as a whole [6] . Each solution 
<&i (r) is determined by the surface charge density on the 
corresponding colloidal species 



(7) 



Given that the cells of each pair of species are consid- 
ered to be neighboring, it is natural to impose a common 
boundary value <&r for the potential on every cell sur- 
face. This value is by definition equal to the average 
value of the electrostatic potential at the cell boundary 
for different species. Hence, we have 



Vt, 



(8) 



V 2 $(r) = 





k 2 sinh$(r) 



for < |r| < a, 
for a < |r| < R, 
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where k — y/ 8ire ou t/ p s (3e 2 is the Debye screening 
length, and $ = j3e$> is the dimensionless electrostatic 
potential, where (3 — l/k^T. This second-order nonlin- 
ear partial differential equation describes the electrostatic 
effects of the ions in the cell volume surrounding the col- 
loidal particle. The boundary conditions are determined 
by the charge on the particle and the cell electroneutral- 
ity The first boundary condition is given by 
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which fixes the electric field at the particle surface. In 
Eq. (H|), Z is the number of elementary charges on the 
particle surface. Because of the homogeneous surface 
charge distribution, there is no electric field inside the 
particle. This reduces boundary condition (f?)) to 



or 



r\.a 
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where Ib = f3e 2 /4ire out is the Bjerrum length. The sec- 
ond boundary condition fixes the electric field at the cell 
boundary, according to Gauss' law 



= 0. 



(6) 



The PB-equation |(5J) together with the boundary condi- 
tions flS} and (JHJ) form a closed set and describe the basic 
Poisson-Boltzmann cell model as studied for example in 
Refs. @, El, p. 

Now we discuss the case of a mixture of equally sized, 
homogeneously charged colloidal species with surface 



where Xi 
that 



Ni/N is the molar fraction of species i, such 



= i. 
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The boundary value <£># is fixed by setting the average 
value of the charge contained in the cells for different 
species to zero. According to Gauss' law, this is achieved 
by imposing 



(10) 



The PB-equation @ applied to each species i together 
with the boundary conditions ([7]), ([8]), and (fTO)) again 
form a closed set. 

In principle, one could allow for a different cell radius 
Ri for each colloidal species. In this factor R? 

must be included in the summation of Eq. (|10p . One 
subsequently imposes a set of physically motivated con- 
ditions on these radii. These conditions must comply 
with the fact that the average cell volume of the system 
is given by 



3 N 



(11) 



In the case of mixtures where the charge of all species 
has the same sign, one can fix each Ri by imposing the 
condition that the electric field must vanish at the cell 
boundary of each species 0]. For reasons of simplicity, 
however, we will not use such an extension. Hence, in 
this paper, the cell radius is given by a single value R. 



III. AN EXTENSION TOWARDS 
HETEROGENEOUS CHARGE DISTRIBUTIONS 

We now consider a system of spherical colloidal parti- 
cles that are not homogeneously charged, such as Janus 
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Figure 2: Illustration of the Poisson-Boltzmann cell model for 
heterogeneous charge distributions. For each direction n per- 
pendicular to the cell surface, we determine the appropriate 
boundary conditions. 

particles Again, the system consists of N identi- 

cal spherical particles of radius a in a volume V. We 
apply the same Poisson-Boltzmann theory, such that the 
electrostatic potential obeys the PB-cquation |(2J) . As be- 
fore, we fix the position of each particle at the center of 
a spherical Wigner-Seitz cell of radius R, and the PB- 
equation must be solved for each "species" . However, in 
this case we replace the species index i by an orientation 
Co, and each solution is denoted by $(u);r). In the spirit 
of Onsager, we can view such an orientational distribu- 
tion as a mixture where the different particle species have 
a distinct charge distribution. 

In this paper, we focus on charge distributions that are 
cylindrically symmetric with respect to the particle orien- 
tation Co. This charge distribution gives us the boundary 
condition on the particle surface at each position an 

$ in (w; an) = $ out (cj; an), (12) 

en-V$m(w;an) = 47W B e r (u);n) + n- V<& out (u>; an), (13) 

where o~(Co; n) is the surface charge density (in units of 
e) that belongs to a particle with an orientation Co, and 
( = fin/tout is the relative dielectric constant of the par- 
ticle with respect to the surrounding solvent. The labels 
"in" and "out" denote the solutions inside the colloidal 
particle, and outside the particle, respectively. 

Next, we must supply a generalized version of the 
boundary conditions on the cell surface given in Eqs. (j8j 
and (flO|) . The following approach is illustrated by Fig. O 
which shows the directionality that must be included in 
the appropriate boundary conditions. First, the concept 
of a fixed cell surface potential is applied to include dif- 
ferent positions on the cell surface. Similar to the cell 
model for mixtures, the value of this potential is defined 



as the average value of the electrostatic potential at the 
cell boundary for different orientations. However, the 
average value is taken at the position of opposite orien- 
tation 

$(u>;Rn) = ($}(-i?n) = $ fl (-n) Vw, (14) 

where instead of a summation over species weighted by 
the molar fractions Xi, we have an integral over particle 
orientations Co weighted by the orientational distribution 
function (ODF) f(u) 

(*)(r) = J dc2>/(w)$(c2>;r). (15) 

This distribution is normalized such that 

J dCof(Co) = 1. (16) 

The definition of the cell surface potential $fl(n) in 
Eq. (|14[) is such that it is always an even function 

* fl (n) = Sfl(-fi). (17) 

Finally, we have to impose a boundary condition that 
fixes this cell surface potential. However, if we only im- 
pose global electroncutrality on the system, we obtain 
a boundary condition that is too general for a solution 
that is not spherically symmetric. It would ensure that 
the average value (over all "species" ) of the charge con- 
tained in each Wigner-Seitz cell vanishes. By applying 
Gauss' law, we see that this condition is satisfied by set- 
ting the average value of the electric field integrated over 
the surface to zero 

J dn[n- V($)(i?n)] =0. (18) 

Interestingly, this fixes only the isotropic contribution to 
the cell surface potential. Therefore, we impose an addi- 
tional condition that is based on the concept of continu- 
ity of the electric field flux from one cell to another. The 
difference between the outward flux at the cell bound- 
ary and the average inward flux of neighboring cells is 
represented by 

AF fl (w; n) = n ■ V$(w; Rh) - n ■ V($)(-i?n). (19) 

This quantity is averaged over all particle orientations, 
and set to zero, in order to insure global electric field 
flux conservation 

(AF fl )(A)=0, (20) 

which is equivalent to imposing 

n V($)(i?n) = n- V($)(-i?n). (21) 

This condition does fix the cell surface potential, and it 
defines an average boundary value of the radial derivative 
such that it is an odd function of n. 



4 



IV. SPECIAL LIMITING CASES: PERFECTLY 
ISOTROPIC AND PERFECTLY ALIGNED 

In this section, we apply specific choices for the ODF. 
In turn, these choices yield a specific form for the bound- 
ary conditions (fl4|) and (j2Tj) . The resulting models are 
less intricate than the full model we presented in the 
previous section. Also, these models yield boundary con- 
ditions that one would expect from a naive description 
of such systems. 

Let us take a look at the model that our boundary con- 
ditions yields when we implement specific ODFs. First, 
we consider a perfectly isotropic orientational distribu- 
tion 

M&) = (22) 

Since in such a system there is no preferential direction, 
we argue that all solutions — for different particle orien- 
tations — are equivalent. Consequently, the cell surface 
potential ^(n) is independent of the position on the 
cell surface 

$ iso (£ ;j Rn) =$ H . (23) 

This result is in accordance with the notion that in the 
isotropic case the average over all particle orientations 
(fT5|) is equal to the average over all orientations n of 
the position on the cell surface, and that this average 
no longer depends on either orientation. Therefore, the 
boundary condition (|2ip is equivalent to the condition 
that each cell is electroneutral 

n V($ iso )(i?n) = fdn[n- V$ iso (w;#n)] = 0. (24) 

Alternatively, we can choose a perfectly aligned orien- 
tational distribution 

f ]] (w)=5(u>-z). (25) 

Clearly, in this case there is only one solution to be de- 
termined 

*ll(i;r) = *||(r), (26) 
and the boundary conditions (|14ll and (f2"T1) read 

$||(i?n) = $||(-i?n), (27) 

n- V$||(i?n) = n- V$||(-i?n). (28) 

Evidently, this choice leads to periodic boundary condi- 
tions. 



V. APPLICATION TO LINEARIZED 
POISSON-BOLTZMANN THEORY 

To solve the full nonlinear problem is possible numeri- 
cally, but it turns out to be very involved [161 ]. Therefore, 



we restrict ourselves to the linearized version of Poisson- 
Boltzmann theory. In this case, the nonlinear right hand 
side of the PB-equation is linearized around a certain 
value. We denote it by $o, such that the linearized 
Poisson-Boltzmann (LPB) equation is given by 

V 2 $ ut(w;r) = K 2 cosh$ ($out(^;r) - $o) + K 2 sinh$ - 

(29) 

In some cases, the value for $o is chosen to be zero. This 
choice is meaningful if the concentration of colloids, as 
well as the total surface charge density, is low. Alterna- 
tively, its value can be set to the isotropically averaged 
value of the potential at the cell boundary. This choice 
is particularly useful when one has the boundary values 
of the potential and the electric field from numerical cal- 
culations of the nonlinear PB-equation [lfjj]. These can 
be used to fit renormalized charge distributions on the 
particle surface using the expression in Eq. ([55)) . Lastly, 
one can apply the Donnan potential as the value around 
which to perform the linearization. This value requires 
no other input than the colloid concentration, its parti- 
cle radius, its total surface charge, and the reservoir salt 
concentration (l4| . In this paper, we leave $o unspeci- 
fied. 

Inside the colloidal particle <E> still satisfies the Laplace 
equation. It is natural in this case to expand both the 
inner and the outer solution in spherical harmonics. This 
leads to two sets of coefficients which have to be matched 
at the particle surface. Inside the particle 

oo +e 

$ in (w;r) =J2 A t ,m{uY tY l,rn{Q,<t>), (30) 

whereas in the cell interior 
$ out (w;r) = $ - tanh$ 

oo +1 

+ Y X! [ B t,m(&)u(Kr) + Ct l7n (ui)kt(Kr)]Yi >m (d,4>), 
e=o m=—l 

(31) 

where R 2 = K 2 cosh$o, and ig and kt are the modified 
spherical Bessel functions of order t of the first and sec- 
ond kind, respectively. The boundary condition on the 
particle surface are given by Eqs. (|12p and (j!3p . where 
we decompose the charge distribution as 

~ 21+ 1 

a(u; n) = ^ er e P e (u; ■ n). (32) 

1=0 

Next, we impose the boundary conditions at the cell sur- 
face, which are given in Eqs. ([T4"[) and ([2"Tj) . Together, this 
yields the general solution for the dimensionless electro- 
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static potential in the cell interior 

$out(w; r) = $o - tanh $ + Z b k V _ _ m 



x [fc^(Kr)i^(Aci?) — w(/€r)fc^(/«iZ)]P^(d) • f) 
__ x ^ (2*+l)<7< 



A^(e; Ka, kR)k 2 R 2 Ei (e; rea, /-ci?) 



£ even 



(33) 



x S^(e; rea, nr) J du' 'f(u')P e (u' ■ f), 
where 

S^(e; Ka, kR) = — k' e (na) — —ki(na) iAkR) 
\ na J 

+ ( iAko) — —ie{na) kAnR), (34) 
\ Ka J 

_ &St(c\ R,a, kR) 



Ai (e; Ka, kR) 



(35) 



d[RR) 

The details of the derivation of the expression in Eq. ([33]) 
can be found in appendix|A] Note that the first sum (over 
all £) does not depend on the ODF, whereas it does de- 
pend on the particle orientation Cj. This contribution to 
the potential is purely due to the particle at the center 
of the cell, and it vanishes at the cell boundary. Con- 
versely, the second sum (over even £) does not depend on 
the particle orientation, whereas it does depend on the 
ODF. This means that it describes the effect of all the 
surrounding particles. Moreover, it vanishes in the limit 
of infinite dilution (i? — > oc). 

The thermodynamic potential for the ion distribution 
in a single cell is given by [3] 

/3ft cc u(w) = p s / dr{$ out (w;r)sinh[$ out (a);r)] 

J out 

-2cosh[$ out (a);r)] +2} 

a 2 f 

+ — I dner(a);n)<I> ou t(w;an), (36) 

where the the label "out" at the integral symbol denotes 
integration over the cell interior (i.e., the domain of $ ut), 
and p s is the reservoir salt concentration (such that k 2 — 
SitIbPs)- We cannot evaluate this expression analytically. 
Therefore, we linearize it around <I>o to find 



where 



47T 



(37) 



/3fi = y(i? 3 -a 3 ) / 9 s ($ sinh$o-2cosh$ + 2), 

(38) 

(3fl iso = p s ($ cosh$ ~ sinh$ ) 

x / dr($ out (^;r)-$ ), (39) 

J out 

a 2 f 

f3tl int = — / dna(Cj;h)<P out (Cj;ah). (40) 



It turns out that /3f2;„t depends on the particle orienta- 
tion and the ODF, whereas the other two terms depend 
on neither. The expression we obtain can be derived 
through another route, by expanding the original non- 
linear functional of the ion profiles p±{Cj; r) up to second 
order with respect to a density /?±,o = Ps exp[=p$o]- Min- 
imizing this functional with respect to the ion profiles 
yields the LPB-equation (f2"9"| , and the accompanying ex- 
pressions for the ion profiles p±(w;r) ~ P±,o[l ±$oT 
<I> out (d;; r)]. Substitution of this expression into the func- 
tional yields Eq. ([37)1 . 



VI. ONSAGER-LIKE SECOND-ORDER 
DENSITY FUNCTIONAL THEORY 

We approximate the total free energy (per colloidal 
particle) of the system by averaging over all particle ori- 
entations. Also, we add an entropic contribution, which 
is analogous to mixing entropy 

j dcj/(w)m[47r/(c2>)] + j dw/(^)/3ficeu(w). 

(41) 

We neglect the translational entropic contributions of the 
colloidal particles, because we are only interested in the 
effects of charge anisotropy and orientational distribu- 
tion. Using the identity 



J dncr(c2;; &)Pi(u)' ■ n) = a e Pi{Cj ■ Cj'), 



(42) 



which can be easily derived from Eq. (|32p using the addi- 
tion theorem, we derive the following expression for the 
free energy 



N 



where 
N 



N 

1 

+ 2 



dw/(«)ln[47r/(o;)] 
dw/(w) / du)'f(u)')K(u>,u)'), (43) 



($ - tanh$ ) 



a 2 On 47T,„, q s 

-^-—(R 3 -a 3 )p s smh$ 



a2<J o ,~ , « \ 

+ ~^-(*o -tanh$ ) 



q 2 Ibk 1 



1=0 



(e; Ka, kR) 



x [kt(Ka)ie(KR) — ii(Ka)ki(KR)] , (44) 
(2£+ l)cr 2 P e (u) ■ Cj') 



and 
K(Cj,Cj') 



E 



£ even 



(45) 



Note that the second term of ^J-q/N vanishes if one 
chooses the Donnan potential for the value of <&o. Also, 
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the entire contribution from 0CFo /N vanishes upon taking 
the functional derivative with respect to f(Co). 

The structure of the free energy functional |43|) is re- 
markably similar to that of Onsager's second- virial the- 
ory for hard rods [l3j. In the case of spherocylinders, the 
kernel K (Co, Co') stems from hard-core interactions, and it 
is equal to the product of the rod density and the ori- 
entation dependent excluded volume of one rod in the 
vicinity of another, 



^-D 3 + 2ttL 2 D + 2i 2 Dsin 7 



, (46) 



where L is the rod length, D is the rod diameter. The 
angle 7 G [0, 7r] between the two rod orientations is de- 
fined by cos 7 = Co ■ Co'. Onsager's second- virial the- 
ory predicts the existence of an isotropic-nematic tran- 
sition, caused by the competition between orientational 
and translational entropy. In the low-density isotropic 
phase, the "mixing" term J dCof(Co) ln[47r/(a>)] is mini- 
mized by an isotropic orientational distribution, whereas 
the contribution due to the average excluded volume 
5 J dCof(Co) J dCo' f (Co') K (Co, Co') is minimized in the high- 
density ncmatic phase. This transition only occurs if 
the lcngth-to-diameter L/D is large, such that the kernel 
K(Co, Co') is sufficiently anisotropic, with a maximum at 
7 = 7r/2. In the limit L 3> D, the description by Onsager 
is quantitative. 

In the present case, the kernel stems from anisotropic 
electrostatic interactions, and we will investigate if these 
can give rise to such a symmetry-breaking transition. In 
both cases, the kernel is rotationally invariant, i.e., it 
only depends on the mutual relative orientation of the 
unit vectors Co and Co' — through the dot product Co ■ Co' . 
Figure [3] shows the values of the kernel K(Co,Co') for two 
distinct surface charge distributions a(Co,h). The angle 
9 between the axis of symmetry and the position vec- 
tor on the surface is defined by cos 9 = Co ■ h. The top 
graph shows a highly peaked distribution around 9 = 
and 9 = n in the inset. However, the kernel is much less 
anisotropic in this case. The lower graph has a purely 
quadrupolar distribution, which reflects in the fact that 
the kernel the same orientational dependence. In both 
cases, the kernel has a minimum at 7 = 7r/2. Conse- 
quently, we expect no isotropic-nematic transition. In 
the next section, we prove that this conclusion holds for 
any choice of parameters. 



VII. BIFURCATION THEORY 

The ODF that minimizes the free energy P5| obeys 
the Euler-Lagrange equation 

1h[4tt/(<&)] = A - J dCo'f(Co')K(Co, Co'), (47) 

where we introduced the Lagrange multiplier A to ensure 
the normalization of / given by Eq. (|16p . The Euler- 
Lagrange equation can be rewritten in a form that always 
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Figure 3: The kernel function K(Co,Co') for different values of 
the cell radius R, and two distinct charge distributions. We 
fixed the values Ra = 1, e = 1, and RIb = 0.01. The inset in 
each graph shows the corresponding charge distribution as a 
function of the angle 9 between the axis of symmetry and the 
position vector on the surface. The charge distribution in the 
top graph is scaled with 00 = 10 3 k 2 , whereas the scaling in 
the lower graph is given by 02 = 10 3 /t 2 . 



satisfies this normalization 

exp[- / dCo'f{Co')K(Co,Co')] 



JdCo' 



exp 



dCo"f(Co")K(Co',Co") 



(48) 



One easily checks that /i SQ (w) = 1/47T is a solution of 
Eq. (|48|) . describing the perfectly isotropic phase. Due 
to its nonlinear character, one can expect additional 
(anisotropic) solutions to this equation. Finding explicit 
expressions for these solutions, however, is difficult, al- 
though good insight can be obtained form a bifurcation 
analysis. The goal of this analysis is to determine if — and 
for what parameters — an instability can be found in the 
reference solution with respect to a perturbation. 

We choose the isotropic ODF ([22]) as a reference, 
and expand around this solution by writing f(Co) = 
1/4-7T + 5f(Co), with Sf a small deviation. Following the 
same scheme as Kayser and Raveche [13] applied to On- 
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sager's model of hard rods — which was extended by Mul- 
der [l8[ — we find the bifurcation equation 



sm = - 




Lu,u')Sf(u') = -IC[Sf}(u). (49) 



This is an eigenvalue equation, for which a non-trivial so- 
lution exists if the integral operator K, has eigenvalues —1. 
The parameter value for which this occurs is called the 
bifurcation point, where an anisotropic solution branches 
off from the (isotropic) reference solution. The solution 
to the bifurcation equation (|49p can be given in terms of 
eigenfunctions of /C. On the basis of rotational-symmetry 
arguments, we find that these eigenfunctions are the Leg- 
endre polynomials of the dot product of the orientation 
to with respect to an arbitrary direction, 

J du>'K(u),u')P e (u}' • z) = \ e P e {uj • z), (50) 

where the eigenvalues Xg follow from Eq. (|4"5)) . 

f 47TCr| 

\ _ J _ r — * —. — ; — - — z — 7=r-, — = — — r-r for i even, 
A( — \ K°a^R' ! Ag(e; na, nR)!Zt(e; na, kR) 

[O for I odd. 

The bifurcation point is determined by Xg = — 1. How- 
ever, all coefficients Xg are positive. Therefore the bi- 
furcation equation (|49p has no solution, and there is no 
bifurcation point. The understanding of the origin of 
this property of the coefficients lies in the fact that both 
S|(e; Ra, RR) and A^(e; Ra, RR) approach their minimum 
in the limit R — ► a. Moreover, these limits are non- 
negative, since 

lim EAe;Ra, RR) = _„ „ , (52) 



lim AAe\Ka, kR) — . (53) 
R^a K A a 6 

This result shows that the breaking of orientational 
symmetry cannot be captured in this simple version of 
Poisson-Boltzmann cell theory. Additionally, it strongly 
suggests that there must be positional order before there 
can be orientational ordering in suspensions of Janus or 
other patchy particles. In other words, the transition 
from an isotropic state to a fully ordered crystal phase — if 
it exists — will be intermitted by a plastic crystal phase. 
Our simple cell model does not take into account the po- 
sitional correlations of the plastic crystal phase, due to 
the mean- field nature of the applied boundary conditions 
at the cell surface. 



VIII. CONCLUSION AND OUTLOOK 

We developed a simple cell model in the context of 
Poisson-Boltzmann theory for heterogeneously charged 



colloidal spheres. The boundary conditions — on the col- 
loid surface as well as on the Wigner-Seitz cell sur- 
face — depend on the charge heterogeneity and the orien- 
tational distribution of the colloidal particles. Within a 
linear approximation to Poisson-Boltzmann theory, these 
boundary conditions give rise to a free-energy functional 
of the orientational distribution function f(uj) that is 
very similar to the one used in Onsager's second-order 
virial approximation in the des crip tion of the isotropic- 
nematic transition of hard rods [13[ . The present descrip- 
tion, however, does not give rise to orientational order- 
ing. Since our model treats the position of the colloids in 
a mean field description — and since we do expect some 
degree of orientational ordering at sufficiently high par- 
ticle density — this result suggests that orientational or- 
dering requires the existence of positional ordering. The 
present theory predicts no orientational ordering in flu- 
ids of these particles, i.e., no liquid crystal phases. On 
the other hand, this effect can also be attributed to an 
oversimplification of our description. In that sense, the 
orientational ordering could be the result of particle pair 
correlations which are not included in our model. 

Therefore, one could consider to expand the class of 
Poisson-Boltzmann cell models even further. Conceiv- 
ably, a method can be devised to include correlations 
of particle positions through the electrostatic boundary 
conditions. More specifically, one can choose a differ- 
ent approach to the way that the surface potential 
is determined. In the present models, this potential is 
independent of the colloidal species to which the cell 
belongs. Also, each colloidal species has an equivalent 
weight — equal to its molar fraction — in the average of 
the potential and electric field flux at the cell boundary. 
This property is due to the mean field description which is 
used. Also inherent to this description is the assumption 
that the surrounding of a particle at the cell boundary 
is independent from the species it belongs to — or equiva- 
lently, its orientation. However, if this restriction is lifted, 
one may include the fact that the surroundings do depend 
on this property through the pair distribution function. 

The nature of these correlations can be related in a 
simple wa y to systems of oppositely charged colloidal par- 
ticles [HI HU- The number of bonds between oppositely 
charged particles in these systems depends on the col- 
loid density. Also, for the dense liquid phase — coexisting 
with a dilute vapor phase provided the Debye screen- 
ing length is large enough — the pair distribution func- 
tion shows that a colloidal particle is surrounded by dif- 
ferent layers of colloidal species with alternating signs of 
charge [19]. The first surrounding layer has an opposite 
charge with respect to the particle in the origin; the fol- 
lowing layer is like-charged, and so on. These systems 
also display multiple crystal structures, which have dif- 
ferent coordination number. The same notion can be sub- 
sequently applied to particles of different orientations to 
include orientational pair correlations in the cell model. 
This paves the way to a description of orientationally or- 
dered phases. Additionally, a jellium approximation can 
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be applied in the same way it is applied to monodisperse 
systems of homogeneously charged colloidal spheres and 
rods [il], I n this description, there is no need for 

a certain cell shape and volume. Moreover, the jellium 
model has a natural way to include particle pair correla- 
tions [23|, Hil • Finally, there is an opportunity to apply 
the Poisson-Boltzmann cell model to non-spherical cells 
[25l | . The boundary conditions can be imposed in the 
same way as in this paper. However, this complicates 
the expression of the appropriate boundary conditions, 
since a non-spherical shape will couple different spheri- 
cal harmonic modes. The shape of these cells must be 
controlled by additional constraints, such as the mini- 
mization of free energy. Also, the choice of shapes must 
be motivated by physical arguments. We leave these op- 
tions for future studies. 



Appendix A: DERIVATION OF THE 
ELECTROSTATIC POTENTIAL 

Inside the colloidal particle $(u>; r) satisfies the Laplace 
equation, whereas in the cell interior it satisfies the LPB- 
equation. Therefore, we have to match two general solu- 
tions, using the boundary conditions on the particle sur- 
face. We do this by expanding both solutions in spherical 
harmonics. These expressions are given in Eqs. ([50)1 and 
(l3"Tj) . We apply the boundary condition on the particle 
surface, given by Eqs. (fT2f and (fT3| . To this end, we 
expand the surface charge distribution in spherical har- 
monics, using Eq. (|32[) and the addition theorem, 



cr(w;n) 



a t YT m (fi)Y itm {n). (Al) 



The arguments n and Cj of the spherical harmonic func- 
tion should be interpreted as the pair of spherical angles 
of this orientation with respect to the reference frame. 
Consequently, from the boundary conditions (TT2"|) and 
([12]) , we obtain the following condition on the coefficients 

Of $out, 

Bf m {u)) (i'ARa) - ^—i e (Ra)] 
\ na J 

+ Ci, m (u>) (k' t (Ra) - ^-ki{Ra)\ 



(A2) 



Next, we apply the boundary conditions at the cell 
surface given in Eqs. (| L4[) and (|2Tj) . This yields a lin- 
ear system of equations, which can be solved analyti- 
cally. However, we can choose to split the solution into 
two contributions. The first contribution then satisfies 
the boundary conditions on the particle surface — given 
by Eq. (|A2|) — as well as the condition that the poten- 
tial vanishes at the cell boundary. This is already the 
relevant boundary condition for all odd contributions to 



5>(<l>;r), whereas a second contribution must be added 
later to the even contributions in order to satisfy the full 
set of boundary conditions. The coefficients that belong 
to the first contribution will be denoted by _B^ m (o>) and 
C^ m {CS). First, we impose the vanishing potential at the 
cell boundary by 

B e , m (tu)i e (RR) + Ci, m {u>)ki(RR) = 0. (A3) 

Together with Eq. (|A2|) . this yields 

Bt, m @) = „ , . f'" h(RR), (A4) 



S^(e; na, kR) 
AttIbR-^iY;^) 
(e; Ra, RR) 



it(nR), 



(A5) 



where S^(e; Ra, RR) is defined in Eq. (|34ll . The orienta- 
tional dependence of this first contribution is such that 
it — and therefore all odd contributions — only depends on 
the angle between cD and f , 

oc 

$ od d(£; r) = IbR- 1 J2 ( 2i + fytPtiw ■ r) 
1=1 

l odd 

ki(Rr)ii(RR) — ig(Rr)kf(RR) 



S^(e; na, kR) 



(A6) 



As previously mentioned, a second contribution must 
be added to the coefficients of the even contributions. 
With this contribution included, the solution $(a>; r) sat- 
isfies the full set of boundary conditions on the cell sur- 
face given in Eqs. (fT4|) and (f2Tj) . We denote the coeffi- 
cients of this secondary contribution by Bi m and Ci >m . 
Also, we show that these do not dependent on the parti- 
cle orientation, because the two distinct boundary condi- 
tions that govern them do not. First, the boundary con- 
dition (|A2[) is already satisfied by the coefficients Be tm (u) 
and Ci >m {ijj). Therefore, 

( 

Bi, m i/(«a) - —ie(Ka) 
\ na 

+ C e . m (k' e {Ra) - J-fa(Ka) ) = 0. (A7) 

\ KCL J 

Second, the boundary condition (|14[) imposes a value on 
the coefficients that only depends on the value of the 
potential at the cell surface, which is necessarily inde- 
pendent of the orientation u>. Hence, 

Bt,mii{RR) + Ce, m ke(RR) = 4>t, m , (A8) 
where (j>g iTn is defined by 

00 +£ 

<P R ; (n) = $ - tanh $ + Yl Yl ^t, m Ye, m (n). (A9) 

e=o m=-t 
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Together, these conditions yield 



B, 



na ) n^(e;fta, kH) 



(A10) 
(All) 



Finally, the boundary condition (|21[) imposes a vanishing 
value of the weighted average of the even contributions 
to the electric field flux at the cell boundary. This condi- 
tion can be expressed in terms of a relation between the 
coefficients B^ m , C^ >m , B( m , and C^ m . By substituting 
the values given in Eqs. (|A4|) and (|A5[) . we arrive at 

B lf7n i' t {RR) + C e , m k' e {RR) 

dw/(w) [B/, m (w)*i(&R) + Ct m {Lu)k' e {RR)] 



Sg(e; na, kR)k 2 R 2 



where we used 



do)/(c2))Y^* TO (cj) for I even, 

(A12) 



fa(RR)i' e {RR) - i l {RR)k' l {RR) = 



R 2 R 2 



W, (A13) 



which can be derived from standard identities for the 
modified spherical Bessel functions. The conditions in 



Eqs. (|A7p and (|A12|) are sufficient to derive similar ex- 
pressions for Bi :in and Cg^ m . However, the construction 
we use to derive the solution to the even contributions en- 
ables us to show that the cell surface potential ^(n) has 
the same symmetry properties as the ODF (in addition 
to the fact that it is composed purely of even contribu- 
tions), 



A^(e; Ra, RR)R 2 R 2 



4a/(a>)r/ ro (w), (A14) 



where Ai(e; na, kR) is defined in Eq. (|35| . With this we 
readily obtain the even contributions 



r)=l B R~ 1 £ 5 



(2£+l)a t 



I even 



S^(e; na, kR) 



x <j [ki(itr)it(KR) - u(Kr)k e (KR)]Pi(u) ■ f) 
S^(e; Ra, Rr) 



Ae(e; na, kR)k 2 R 2 



dw7(w')P/(w' • f) , 



(A15) 



and we obtain the general solution for the dimension- 
less electrostatic potential in the cell interior given in 
Eq. ~ 
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